Use case: Tracking water level with GeoMet-OGC-API¶

# Module importation
from datetime import date

from dateutil.relativedelta import relativedelta
import pandas as pd
from owslib.ogcapi.features import Features
from osgeo import ogr, osr
import plotly.offline as pyo
pyo.init_notebook_mode()

# Parameters

# Coordinates
lat = 49.162676
long = -121.958943

# Buffer size in kilometer
buffer = 100

# Number of months used to analyse historical data
numMonths = 3

# Prefered projection
# NAD83 / Statistics Canada Lambert
projection = 3347

# Generating parameters for the OGC API - Features request

# Bounding box a little bigger than buffer size (1 km ≈ 0.01°)
bbox = [str(long-(buffer/2)*0.01),
        str(lat-(buffer/2)*0.01),
        str(long+(buffer/2)*0.01),
        str(lat+(buffer/2)*0.01)
       ]

# Time
endDate = date.today()
endDate = date(2013, 7, 20) # Available data to test the code, will not be in final use case
startDate = endDate - relativedelta(months=numMonths)
startDate = startDate.strftime("%Y-%m-%d")
time = f'{startDate}/{endDate}'

# Retrieval of hydrometric stations data
oafeat = Features('https://geo.weather.gc.ca/geomet/features')
stationData = oafeat.collection_items('hydrometric-stations', 
                                      bbox=bbox, 
                                      STATUS_EN='Active')

# Generating a list of stations inside the buffer zone

# Set the driver
driver = ogr.GetDriverByName('GeoJSON')

# Create a datasource and get the layer
stationData = str(stationData).replace("'", '"').replace("None", "null")
dataSource = driver.Open(stationData, 0) # 0 means read-only. 1 means writeable.
layer = dataSource.GetLayer()

# Identify input spatial reference system
srsInput = layer.GetSpatialRef()
sr = osr.SpatialReference(str(srsInput))
epsg = sr.GetAuthorityCode(None)
srsInput.ImportFromEPSG(int(epsg))

# Spatial reference system used to project data
srsProjected = osr.SpatialReference()
srsProjected.ImportFromEPSG(projection)

# Transformation from input srs to NAD83 / Statistics Canada Lambert
transform = osr.CoordinateTransformation(srsInput, srsProjected)

# Create a buffer to select stations
point = ogr.Geometry(ogr.wkbPoint)
point.AddPoint(long, lat)
point.Transform(transform)
pointBuffer = point.Buffer(buffer*1000) # The value must be in meters

# Select the stations in the buffer zone
stations = []
for feature in layer:
    geom = feature.GetGeometryRef().Clone()
    geom.Transform(transform)
    if geom.Intersects(pointBuffer):
        stations.append(feature.GetField('STATION_NUMBER'))

# Temporary prints, will not be in the final use case
print(len(stations))
print(stations)

# Retrieval of hydrometric data for each station

# Dictionary that will contain a data frame for each station 
waterLevel = {}

# Retrieval of the data and creation of the data frames
for station in stations:
    hydroData = oafeat.collection_items(
        'hydrometric-daily-mean',
        bbox=bbox,
        time=time,
        STATION_NUMBER=station
    )
    dates = []
    levelValues = []
    for feature in hydroData['features']:
        dates.append(feature['properties']['DATE'])
        levelValues.append(feature['properties']['LEVEL'])
    waterLevelData = {'Daily mean': levelValues}
    dfWaterLevel = pd.DataFrame(waterLevelData, index=dates, columns=['Daily mean'])
    waterLevel[station] = dfWaterLevel

# Creation of interactive plots

pd.options.plotting.backend = 'plotly'

# Function to add constant functions to the plot
def constantLine(station, yValue, varName, lineColor='black', lineType='solid'):
    fig.add_scatter(y=[yValue] * len(waterLevel[station].index),
                x=waterLevel[station].index, 
                mode='lines',
                name=varName,
                line=dict(color=lineColor, 
                          dash=lineType)
                )
    
    
station = '08MH147'

# Creation of the plot
fig = waterLevel[station].plot(title='Water level at station {}'.format(station), 
                               template='plotly_white',
                               labels=dict(index='Date', 
                                           value='Water level (m)', 
                                           variable='Variable'),
                               color_discrete_sequence=['gray'])

# Adding the historical maximum and minimum and the current value
constantLine(station, 
             float(waterLevel[station].max()), 
             'Historical max', 
             'red', 
             'dash')
constantLine(station, 
             float(waterLevel[station].min()), 
             'Historical min', 
             'limegreen', 
             'dash')
constantLine(station, 
             waterLevel[station].iloc[-1][0], 
             'Current value', 
             'cornflowerblue')

# Updating some plot parameters
fig.update_traces(hovertemplate=None)
fig.update_layout(hovermode='x unified')
fig.update_layout(yaxis=dict(range=[float(waterLevel[station].min() - 0.1),
                                    float(waterLevel[station].max() + 0.1)]))
# fig.update_xaxes(tickformat="%Y-%m-%d")

pyo.iplot(fig, filename = 'offline')